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1. Introduction 

Non-perturbative calculations of the phases strongly interacting matter are hampered by the 
sign problem (for a review of different approaches see [||]). 

In the past years, however, progress in complex Langevin (CL) dynamics ^ ^ has led 
to the hope that the full region of physical interest can be explored. In this study we support the 
complex Langevin simulations by exploring an alternative approach to lattice fermions at nonzero 
chemical potential: the hopping parameter expansion, which can be formulated as a systematic 
approximation to finite density QCD. The expansion is expected to converge at not too small quark 
masses. 

Historically the hopping expansion was used in the form of the loop expansion (described in 
detail in Sec. 3). It was used by earlier studies at LO and NLO ^ with the full Yang- 

Mills action, also to map the phase diagram in [|^, |^. In studies combined with strong coupling 
expansion it has been possible to calculate NNLO contributions as well 1^ [TJ, but it has 
proven quite difficult to further extend the the expansion to higher orders. 

Here we present an alternative way to introduce higher-order corrections in the hopping pa¬ 
rameter expansion The approach allows calculations at very high orders (only limited by 

available computer power), while keeping the full Yang-Mills action, and without having to con¬ 
sider fermionic loops and their combinatorial factors at each new order. 

We define fhe K- and fc^-expansions below, wifh slighfly differenf properties. We improve on 
fhe convergence properfies of fhe loop expansion where fhe effective expansion paramefer is kNt: 
wifh Nt fhe femporal exfenf of fhe laffice. 

In section 2, we briefly describe Complex Langevin simulafions. In Section 3, we firsf review 
fhe loop expansion, fhen describe fhe new approaches we call K- and tCs-expansion, and discuss 
fheir implemenfafion in fhe complex Langevin equafion. In Section 4 we presenf numerical resulfs 
gained using fhis approach. Finally, in Secfion 5 we conclude. 


2. Complex Langevin Simulations 


The Complex Langevin approach is based on fhe complexificafion of fhe Langevin equafion 
10 0- This also leads fo fhe complexificafion of fhe field manifold. The resulfing process is 


susceptible fo numerical problems (runaway frajecfories, solved by using adaptive sfep sizes [18]), 
as well as convergence fo a wrong resulf. Recenfly if has been shown fhaf convergence is guaranfeed 
as long as some condifions are safisfied, such as fhe fasf decay of field disfribufions and holomorphy 
of fhe action and fhe observables [^. Nofe fhaf fhere are several types of modifications possible 
fo adapf fhe Langevin process for a given action, which one can use fo gef fhe process fo satisfy 
convergence criferia [Hi- The mefhod has proven useful in ofher contexfs wifh a complex action 
as well 


In laffice QCD fhe discrefised Langevin equafion is wriffen as [125|] 


^ exp Y,^Xa{eKj,ya + Ver] X 


a 


T.V) 


(2.1) 
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where Kxva = —D^vaS is the drift force, e the (adaptive) stepsize, and rj independent Gaussian 
noises satisfying (rjxva’n.x'v'a'} = '^Saa'Sxx'Sw'■ A complex action leads to a complex drift K, and 
link variables take values in SL(3,C), losing their unitarity. 

The available configuration space is thus complexified, and loses compactness. For gauge 
theories this leads to an additional complication: the volume of gauge orbits corresponding to a 
configuration is infinite. To restrict the movement of the system along the infinite gauge orbits 
one has to modify the process, while respecting the gauge invariance of the action and observables. 
This can be very conveniently achieved with the gauge cooling (see also the review [^]), which 
uses non-compact gauge transformations to force the process to stay near the unitary manifold. As 
a consequence the decay of the distributions is fast, as required for the convergence proof. Together 
with the adaptive stepsize this practically eliminates runaways. 

Another requirement for the proof of convergence is the holomorphy of the action, which is 
unfortunately not satisfied for QCD. This manifests in zeros of the measure, i.e. detM = 0, leading 
to a meromorphic drift. Poles in the drift then might lead to wrong convergence of the process, as 
shown in nontrivial, soluble models |^], while in many cases the process gives correct results in 
spite of a non-holomorphic action, especially in the cases where the non-holomorphy is due to a 
Haar measure or Jacobian [Hi- 


3. Hopping parameter expansions 

3.1 Loop expansion 

Recall the path integral formulation of QCD, where we use the plaquette action S'ym for the 
gauge fields 

Z = jDUt^^, 5 = 5YM-logdelM, (3.1) 

with the Wilson fermion matrix M, the hopping term Q of which we split into spatial hopping terms 
S and temporal hopping terms R 


M = 1-kQ=1-k,S-R, (3.2) 

with 

Sxy = 2Y^ +r+iU-,^dy,x-a ), (3.3) 

i=l ^ ' 

Rxy = , 

using the matrices r±v = (1 =t 7v)/2. Note that these matrices are projectors satisfying F^^, = r±v 
and r+vL-v = 0. We then expand the fermionic part of the measure 

ff" 

—Tr ^ — Q" 

n=\ ” 

which we can rewrite noticing that we can perform the sum for each loop built from hopping terms 
on the lattice separately, then resum the determinant for each loop (ignoring possible convergence 
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T =1/N^a, 


Figure 1: An illustration of the different loops contributing to the loop expansion of the fermion determinant. 


problems) 

—Tr E —^c =ndet(l-»cCL,), (3.5) 

C,s=l ^ ) c 

where C goes over all the possible loops on the lattice and Ic is the length of the loop. Since 
r+yP-v = 0, no loop that turns back on itself needs to be considered. 

In the static limit, where tc —)■ 0, —)• oo, (^ = 2k: exp= const., all contributions vanish except 
for the Polyakov loops. This is called the leading order of the loop expansion. Next to leading 
order (NLO) loops are gained by a decorating the Polyakov loop with two spatial hoppings, as 
illustrated in Fig. ||. NNLO contributions involve Polyakov loops with more decorations as well as 
the plaquette contribution. 

3.2 K expansion and fc, expansion 

To define an expansion which can be conveniently included in the Langevin equation, we go 
back to the first step in the loop expansion 

( °° k” \ 

—t^E —j • 

In the sum only even powers of the hopping matrix contribute, as a trace is present. The expression 
is straightforwardly generalized to Nf > 1 flavors. We call fhis expansion fhe tc-expansion. The 
fermionic observables can fhan be expressed in terms of fhe expansion as 

iw) = (3-7) 

for fhe chiral condensate (yy) baryonic densify n, using = N^N-i fhe lattice volume. 
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Alternatively, we can rewrite the fermion matrix using the identity 

M = {\-R)(^l--^K,Sy (3.9) 

before expanding to gain the following expansion: 

“ k " / I \" 

detM = det(l —/?)exp ^-y l— R^J ’ (3.10) 

Since the matrix (1 —7?)^^ is diagonal in spatial indices (and dense in temporal indices) and S 
describes spatial hoppings, again only even powers contribute in the sum. The first factor of the 
expansion det(l —R) describes the LO contribution in the static limit, also known as the HDQCD 
approximation. This simplification happens only for Wilson fermions, as backtracking is forbidden, 
hence the only possible loop built from temporal hoppings is the Polyakov loop. In this case the 
determinant and the inverse of the matrix can be calculated analytically as follows. The inverse of 
the matrix can be written as 


(1 -R)-^ = {\-R+y^ + {\-R-y^ -1 ( 3 . 11 ) 

with R^+R = R, R^R = 0 and /?+ (R^) containing hoppings in positive (negative) temporal 
directions. We can then find the inverse of the two terms by the simple expansion (omitting spatial 
coordinates) 


(1 -/?+)^i = £ (2fcc^^r_4f4,45v„v+«4)", (3-12) 

n—0 

Separating the parallel transporter between x and y, we can easily resum the remaining factor to 
give 


;i-/?+)^‘ = i-r_4 


{lKe^Y^P{x) 
l + {lKei^Y-P{x) 


ifv = y 


(3.13) 


= —r_4(2K'C 


l + {2Kei^Y^P{x) 

1 


g\NT:+y-x 


1 + (2k-cA')^^P(v 


W{x,y) ify>x 

-W{x,y) ify<x 


where W(x,y) is the parallel transporter between v and y built from positive hoppings, P{x) is 
the untraced Polyakov loop starting from site x (that is P{x) = W{x,x)). The inverse of (1 —7?^) is 
calculated similarly. Similar formulas were derived in Refs. [ |T^ ] in an effort to develop an effective 
theory for Polyakov-loops, also utilizing the strong coupling expansion for the Yang-Mills action. 

The observables in the fQ expansion can be recovered using the defining equations such as 
{n) = dlnZldpi in sfraighf for ward calculations fo yield formulas similar fo (3.7). 

These fwo expansions can be very convenienfly implemenfed in Langevin simulafions, as de- 
failed in fhe nexf subsecfion, buf fhey have differenl sfrengfhs and weaknesses. We consider a 
fruncafed version of fhe expansion fo order N^LO, in which ferms up fo are kepf. We keep also 
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the terms proportional to exp(—/i), although their contribution is suppressed at large ju,, but they 
lead to the determinant satisfying the symmetry: 

detM(At) = {det{M{-n*))* (3.14) 

The K is expansion is very cheap to calculate, but its convergence properties at nonzero chemi¬ 
cal potential are not optimal, as Q has terms proportional to K:exp(/i). The action truncated to some 
order is holomorphic, so proofs of convergence (requiring also fast decaying distributions) apply 
1^]. In the K expansion one needs to go to order in order to see any K dependence, as the first 
closed loop (the plaquette) can be formed at this order. Similarly, one needs to go to order in 
order to see /r dependence (using the Polyakov loop), since for shorter loops the dependence can¬ 
cels. As we will demonstrate below calculating high orders is easy in this setup, so this drawback 
is not a serious one. 

In our second scheme, the tc^-expansion the main part of the /r dependence is dealt with ana¬ 
lytically, therefore one expects better convergence properties at high /r, and this is indeed satisfied, 
see in Section. The price to pay is the slightly higher numerical cost and the non-holomorphic 
action. 


3.3 K- and expansion in Langevin simulations 

It is useful to consider the expanded determinant ( ^^ as part of the action. Note that the 
resulting effective action is holomorphic. It has than a contribution to the drift term of the Langevin 
equation ( |2.1| ) 

oo 

K,ya = - £ tc'-Tr {Q’'-^D,yaQ) . (3.15) 

n=l 

Note that this contribution is non-real, therefore we have to complexify the theory and use complex 
Langevin dynamics. (At /r = 0 real Langevin simulations are possible by taking the real part of the 
fermionic drift terms.) 

This contribution to the drift term can be estimated using a random vector T], (where / repre¬ 
sents space-time, colour and Dirac indices), with the properties (t],) = 0, (t]f t]/) = as 

K,ya = {V*{D,yaQ)s), ^ T]. (3.16) 


The calculation of this term at N^LO thus requires 2q multiplications with the sparse matrix Q in 
every timestep, when the random vector is refreshed. In the case of the tQ-expansion, the drift term 
has contributions from several places. The lndet(l —R) term has contributions to the drift identical 
to the HDQCD, as calculated in |28]. The contribution of the expansion is 


K,ia = - ^ <Tr 


1 


n—\ 


K,Aa = -Y. <Tr I 


n—\ 


l-R 

Th 


(D^iaS) 

{DxAaR) 


1 


-| n—\^ 


l-R 

1 

1 -/?' 


(3.17) 


for spatial and temporal links, correspondingly. These contributions are estimated using noise 
vectors similarly to (3.16). The numerical effort of the expansion involves also multiplications 
with (1 — 7?)^\ calculated according to (3.13), as well as multiplications with the sparse S matrix. 
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Figure 2: Dependence of the quark density (in lattice units) on the order of the truncation of the K- and 
tfv-expansions, for ji = 0.7 and 1.1, on a 4^* lattice with j3 = 5.9, K = 0.12, and Nf = 2. The region where 
the tc-expansion breaks down for /i = 1.1 is indicated. The lines show the result for full QCD. Saturation 
density is Wsat = ^NcNf = 12. 
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4. Numerical results 


We have used simulations on small 4^ lattices to examine the convergence properties of the 
expansions. We have used two flavors of Wilson fermions with jS = 5.9 and K = = 0.12, for 

several values. We compare results with full QCD, obtained by complex Langevin simulations, 
extending the previous results available for staggered fermions [Q] to Wilson fermions. Since the 
full QCD result is also obtained with Wilson fermions there is no need to renormalize the results, 
already the bare quantities of the expansions in lattice units should converge to their full QCD 
values. We plot results in lattice units. The lattice spacing is measured using the gradient flow, as 


proposed in Ref. [[29|]. The lattice spacing depends weakly on the mobility of the fermions for the 
heavy quark masses that we are using. For HDQCD we find that j3 = 5.9 and K = 0.12 corresponds 
to a ~ 0.12 fm, while for full QCD we find a ~ 0.114 fm. 

In Fig. 1^ we show the quark number density as a function of the order of the K- and fQ- 
expansions for two different /r values, comparing to the full QCD result. We see that at the smaller 
chemical potential /r = 0.7 both expansions behave similarly, with convergence around the 
order. At the higher chemical potential value, the K expansion breaks down, as expected. 

We observe good convergence of the series to the full QCD values. This is a non-trivial 
agreement which supports both the expansion and the full QCD. This means in particular in the 
case of the K expansion, that we must obtain the correct value also in full QCD. This apparently 
means that the non-holomorphicity of the action for full QCD is not a problem (at least for the case 
where the K expansion converges). 

In Fig. 1^ we see similar behavior for the chiral condensate and the spatial plaquette average. 
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Figure 3: Dependence of the chiral condensate and the spatial plaquette average (in lattice units) on the order 
of the truncation of the K- and (Cj-expansions, for p = 0.7 and 1.1, on a 4^ lattice with = 5.9, K = 0.12, 
and Nf = 2. The region where the ff-expansion breaks down for /i = 1.1 is indicated. The lines show the 
result for full QCD. 
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Figure 4: The density at if = 0.14 on a 4“* lattice for the ifs-expansion and full QCD. 


While the convergence appears quite quick at K = 0.12, at larger K this might not be the case, 
see Fig. ^ where we show the performance of the Kg expansion at fc = 0.14. These results suggest 
that the convergence radius of the Kg expansion is below fC = 0.14 at /x = 0.9. The convergence 
radius seems to be independent of the lattice size, however. In Fig. we show the convergence of 
the density on a 8^ lattice. This lattice system has a temperature below the deconfinement transition. 
We see that density grows about a factor of 3 as one changes the chemical potential from p = 0.7 
to p = 0.8, which is a sign of the rapid onset transition. One observes that the Kg expansion still 
performs well in this region. 

5. Conclusions 

In this study we have presented the novel implementations of the hopping expansion for the 
complex Langevin equation which are called K- and Kg expansion. They allow calculations at 
very high, previously impossible orders. This allows the direct observation of the convergence of 
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Figure 5: The density at ic = 0.12 on a 8^ lattice for the JQ-expansion and full QCD for chemical potential 
fi = 0.7 and fi = 0.8. 


the series to the full QCD result. We use the complex Langevin equation to circumvent the sign 
problem of these theories at finite chemical potential. We use no further approximation other than 
the hopping expansion, the gauge action is kept intact in particular, thus our method can also be 
used as a test ground for other effective models. 

Our expansions have different merits: the K expansion is cheap and has a holomorphic action, 
but its convergence properties are bad at large chemical potentials. The Kg expansion is slightly 
more expensive numerically, but has improved convergence properties also at high chemical poten¬ 
tials. 

We performed simulations of the expansions and observed good convergence to full QCD at 
not too high K parameters. This convergence supports both the expanded and the full theory, as 
the agreement shows that the non-holomorphy of the full theory has apparently no impact on the 
results, at least in the parameter range where the convergence is observed. 

The first results indicate that at least the onset transition might be within the reach of this 
method in the cold and dense region of the QCD phase diagram, but further studies are required at 
low temperatures on large lattices. 
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